Zone Segmentation V2

This model was last run on 2020-04-01 16:11:40

Contents

Background and Scope

Feature Selection

Clustering - Kmeans

Bootstrapping - Kmeans

Final Outputs

Appendix

Clustering - Hierarchical Clustering

Bootstrapping - Hierarchical Clustering

Background and Scope

Different teams at Deliveroo have come up with different zone segmentation methods in the past. However, all of these have been made with a specific agenda. For example, Project Hopper has segmented zones based on the costs required to improve their selection-service trade off, while the Pricing team has segmented zones based on price-related factors.

The aim here is to come up with a non-strategic set of zone clusters, which clusters zones according to their fundamental features.

The scope of this analysis is to segment zones in RUK, before automating the segmentation method across different countries. All features were computed using data between 1st January and 29th February 2020. This was to avoid abnormal behaviour due to the Coronavirus outbreak. Competitor features were computed using data between 1st and 7th March 2020, due to unreliable Ubereats scrapes before March.

#df$comp_uber_rx_supply[is.na(df$comp_uber_rx_supply)] <- 0
df$pct_mplus <- df$rx_mplus/df$live_rx_total
df$aof <- df$aof/2 #silly sql mistake sorry

df$zone_tenure_days[df$zone_code == 'LVN'] <-  as.numeric(as.Date('2020-03-14') - as.Date('2020-01-09'))

#df$comp_number_competitors <- df$comp_justeat_bool + df$comp_uber_bool


factors_demo <- c('zone_tenure_days', 'zone_population', 'area', 'popn_density', 'zones_in_city', 'pct_zone_of_city')
factors_rx_supply <- c('live_rx_total', 'rx_per_thousand_popn', 'has_editions_site', "rx_core", "rx_editions", "rx_mplus", "pct_mplus", 'avg_price_category', 'number_cuisines')
factors_rx_cuisine <- names(df)[grepl("num_.+_rx", names(df))]
factors_rx_type <- names(df)[grepl("num_rx.+", names(df))]
factors_comp <-  names(df)[grepl("comp.+", names(df))]
factors_cust <- c("orders", "customers", "customers_new", "customers_existing", "aof", "aov", "avg_credits_vouchers", "plus_orders_pct", "convenience_orders_pct")
factors_riders <- c("cnt_riders_3w", "scooter_deliv_pct", "bicycle_deliv_pct", "car_deliv_pct")

live_rx_cc <- rowSums(df[,factors_rx_cuisine])
for(i in factors_rx_cuisine) {
  df[[paste0(i,"_pct")]] <- df[[i]]/ live_rx_cc
}

live_rx_type <- rowSums(df[,factors_rx_type]) 
for(i in factors_rx_type) {
  df[[paste0(i,"_pct")]] <- df[[i]]/ live_rx_type
}

df$num_rx_enterprise_pct <- df$num_rx_entgold_pct + df$num_rx_entsilv_pct
df$num_rx_local_pct <- df$num_rx_logold_pct + df$num_rx_losilv_pct + df$num_rx_lobron_pct

factors_rx_cuisine <- names(df)[grepl("num_.+_rx_pct", names(df))]
factors_rx_type <- names(df)[grepl("num_rx_.+_pct", names(df))]

Feature Selection

We group features into the following segments:

1. Demographic features: Zone area, population density, number of zones in city, etc
2. Restaurant supply features: Number of live restaurants, number of core restaurants, zone average price category, etc.
3. Restaurant cuisine features: Number of restaurants available for each Deliveroo cuisine category
4. Competitor features: Booleans for which competitors are active in the zone, total restaurant supply per competitor
5. Customer features: New customers, existing customers, AOF, AOV, percentage of plus orders, etc.
6. Rider features: Number of total riders, percentage split by vehicle, etc.

We do not use features like EOD and number of B10s, as these are levers that can be easily pulled by Deliveroo, rather than underlying features of a zone. We ultimately narrowed these features down into 12 variables that we used, to leep the model parsimonious.

We normalize each variable, and use correlation matrices to remove features that are highly correlated to each other. Below is the correlation matrix for restaurant-related factors.

df_vars <- df[,c(factors_demo, factors_rx_supply)]


normalize <- function(x) {
return ((x - min(x)) / (max(x) - min(x)))                        
}


cols_to_normalize <- names(df_vars)[1:ncol(df_vars)]

for (i in cols_to_normalize) {
  df_vars[[paste0("norm_", i)]] <- normalize(df_vars[[i]])
}

df_vars2 <- df_vars[,-which(names(df_vars) %in% cols_to_normalize)]
names(df_vars2) <- gsub("norm_", "", names(df_vars2))


c <- cor(df_vars2[,1:ncol(df_vars2)])

corrplot(c, method = "number", type="upper", order="original",
         col=brewer.pal(n=8, name="RdYlBu"), tl.cex =0.8, number.cex = 0.6)

NA
NA
NA

Here is a correlation matrix with the main demographic and customer features.

We thus remove features with a correlation score above 0.7. These are: percentage population of zone vs city , total core restaurants, number of mplus restaurants, existing customers, etc. For restaurant type, we only keep features for G6, combined enterprise (Gold, Silver) and combined local restaurants.

Just a note, for variables around rider vehicle type, some zones have all vehicles classified as NULL. These are predominantly Marketplace-only zones.

df_vars <- df[,c("zone_code", factors_demo, factors_rx_supply, "num_rx_g25_pct", "num_rx_enterprise_pct", factors_comp, factors_riders)]

cols_to_remove <- c('pct_zone_of_city', 'rx_core', 'orders', 'customers', 'customers_existing', 'rx_mplus', 'rx_editions', 'num_african_rx_pct', 'num_malay_rx_pct', 'cnt_riders_3w', 'rx_per_thousand_popn', 'avg_price_category', 'scooter_deliv_pct', 'bicycle_deliv_pct', "num_rx_local_pct", "has_editions_site", 'comp_justeat_bool',
                    'comp_uber_bool', "area")
df_vars <- df_vars[,-which(names(df_vars) %in% cols_to_remove)]
        
cols_to_normalize <- names(df_vars)[2:ncol(df_vars)]

for (i in cols_to_normalize) {
  df_vars[[paste0("norm_", i)]] <- normalize(df_vars[[i]])
}

df_vars2 <- df_vars[,-which(names(df_vars) %in% cols_to_normalize)]
names(df_vars2) <- gsub("norm_", "", names(df_vars2))
#df_vars2 <- subset(df_vars2, select = - c(num_african_rx_pct, num_malay_rx_pct))

We do a quick Principal Component Analysis (PCA) to understand which features are correlated, and also see what differentiates different zones from each other. The first 2 principal components explain 0.60309 of the total variation, so do keep in mind that there are other orthogonal principal components we have not visualised.

Below is a biplot based on our dataset.

set.seed(1234)
#df_vars2 <- df_vars2[1:10,]
df.pca <- prcomp(df_vars2[,2:ncol(df_vars2)], center = T, scale = T)
comp <- cbind(zone_code = df_vars2$zone_code, data.frame(df.pca$x))


rownames(df_vars2) <- df_vars2$zone_code
#summary(df.pca)

#biplot(df.pca)

# autoplot(df.pca, shape = F,
#         data = df_vars2, label = T, label.size = 2,
#         loadings = T, loadings.label = T, loadings.colour = 'blue', loadings.label.size = 4,
#         loadings.label.vjust = 0, loadings.label.hjust = 0,
#          xlim = c(-0.06,0.165), ylim = c(-0.15,0.175), position = "jitter")

autoplot(df.pca, data = df_vars2, label = T, shape = F, label.size = 2.5,
         loadings = T, loadings.label = T, loadings.colour = 'blue', loadings.label.size = 4,
loadings.label.vjust = -1, loadings.label.hjust = 1)

NA
NA

Factors that are most important for differentiating zones are number of live restaurants in the zone, number of cuisines, how long the zone has been open (zone_tenure), the percaentage of M+ restaurants, JustEat and Uber’s restaurant supply, etc.

Some factors are significantly correlated to zone size and density, such as UberEat’s restaurant supply.

The biplot allows us to see what makes 2 zones similar or dissimilar. For example, MC (Manchester Central) and BNC (Brighton Central) are pretty similar along features that affect PC1 - high zone tenure, high number of cuisines, high restaurant supply, both have Editions sites, etc. What differentiates them is the number of zones in their respective city/metropolitan areas (26 vs. 1 zone), percentage of car deliveries in the zone (1% vs. 11%), their competitive landscape (Justeat and Ubereats have much higher supply in MC, whereas Deliveroo supply is roughly equal in both zones), etc.

Clustering using Kmeans

We now cluster this data using kmeans. A scree plot shows that the optimal number of clusters is around 5 or 6, as beyond that, the reduction in the total within-cluster sum of squares tapers off.

k.max <- 15

wss <- sapply(1:k.max,
              function(k){kmeans(df_vars2[,2:ncol(df_vars2)], k, nstart=50,iter.max = 15 )$tot.withinss})
plot(1:k.max, wss,
     type="b", pch = 19, frame = FALSE,
     xlab="Number of clusters K",
     ylab="Total within-clusters sum of squares")

We then segment the zones into 5 clusters, seen below.

seed <- sample(c(1:1000),1)
set.seed(926)


k <- 5
rownames(df_vars2) <- df_vars2$zone_code
op.kmeans <- kmeans(df_vars2[,2:ncol(df_vars2)], k)


clusters <- data.frame(zone_code = rownames(as.data.frame(op.kmeans$cluster)),
                       cluster = op.kmeans$cluster)


df_final <- merge(clusters, df,
                  by = 'zone_code',
                  all.x = T)


autoplot(op.kmeans, data = df_vars2, label = TRUE, label.size = 3
         , shape = F
         , frame = T)

#print(seed)

The clusters represent 5 groups of zones: Urban Winners, Outer City, Momentum - New Launch and Relaunched, Opportunity - Wealthy and University Towns, and Why are we here? - M+ Zones. Below is a summary of their main characteristics.

clus1 <- c('Nickname: Urban Winners'
, 'Very large zone population and very high population density'
, 'Zone is >4 years old'
, 'Mix of large and smaller metropolitan areas - mainly city centres'
, 'Highest rx supply, very high cuisine diversity, not many M+ restaurants'
, 'Deliveroo dominates the competition'
, 'Only cluster with Editions sites, not many M+ restaurants'
, 'Only cluster where most deliveries are on bikes/scooters'
, 'Example zones: Brighton Central, Birmingham, Cambridge, Liverpool City Centre, Edinburgh North/South, Portsmouth'
, '')


clus2 <- c('Nickname: Opportunity'
, 'Medium zone population'
, 'Most zones are ~ 3 years old'
, 'Small metropolitan areas, most zones are the only zone in their city'
, 'Decent rx supply, low number of M+ restaurants'
, 'Good proportion of G6 and Enterprise restaurants'
, 'Relatively high density of restaurants, good cuisine diversity'
, 'Justeat and Deliveroo competing, Ubereats present but has lower rx supply'
, 'Example zones:  Aberdeen, Cheltenham, Slough, Lancaster, Exeter, Warwick, Durham'
, '')


clus5 <- c('Nickname: Outer City'
, 'Large zone population'
, 'Zone is 2 - 3 years old'
, 'Very large metropolitan areas - most zones are outer Birmingham and Manchester'
, 'Decent rx supply, decent cuisine diversity'
, 'Justeat dominates in terms of restaurant supply'
, '37% (second highest) M+ restaurants'
, 'Example zones: Walsall,  Bolton, Wolverhampton, Dudley, Edgbaston'
, ''
, '')


clus4 <- c('Nickname: Momentum - New Launch and Relaunched'
, 'Medium zone population and low population density'
, 'Most zones < 1 year old, some older zones'
, 'Mainly new launch or relaunched (eg Livingston) zones, handful of older zones (Bangor, Beeston, Bedford)'
, 'Low rx supply, low cuisine diversity'
, 'Justeat doing better than Deliveroo and Uber - Uber not in 25% of zones.'
, 'Highest proportion of G6 restaurants, very high M+ restaurants'
, 'Mainly car deliveries'
, 'Example zones:  Doncaster, Scunthorpe, Weymouth, Sunderland, Burton upon Trent, Derry City Centre'
, '')


clus3 <- c('Nickname: Why are we here? and M+ only zones'
, 'Range of zone populations, and medium population density'
, 'Range of small to large metropolitan areas'
, 'Most zones are < 2 years old'
, 'Lowest rx supply, lowest cuisine diversity'
, '84% M+ restaurants'
, 'Justeat dominates, Uber not present in 20% of the zones'
, 'Main difference between this cluster and << is M+ restaurants, population density, and justeat rx supply'
, 'Example zones: Outer Newcastle, Outer Rotherham, Outer Sheffield, Manchester East, Seaham'
, '')



clus_desc <- data.frame(Cluster1 = clus1,
                   Cluster2 = clus2,
                   Cluster3 = clus3,
                   Cluster4 = clus4,
                   Cluster5 = clus5)

clus_desc

NA

Here is a more detailed view of the features in each zone. Take a look at the biplot above to understand which features are more important in differentiating them. The most important features are number of live restaurants, number of cuisines, zone tenure, competitive supply and percentage of M+ restaurants. We’ve also summarised some features that were not used in the clustering, just to illustrate the difference in the clusters more clearly.

df_summ <- df_final %>%
  group_by(cluster) %>%
  summarise(number_of_zones = n(),
            zone_tenure_days = round(mean(zone_tenure_days)),
            zone_population = round(mean(zone_population)),
            #area = mean(area),
            popn_density = round(mean(popn_density)),
            zones_in_city = round(mean(zones_in_city)),
            live_rx_total = round(mean(live_rx_total)),
            number_cuisines = mean(number_cuisines),
            #rx_per_thousand_popn = mean(rx_per_thousand_popn),
            has_editions_site = mean(has_editions_site),
            pct_mplus = mean(pct_mplus),
            avg_price_category = mean(avg_price_category),
            car_deliv_pct = mean(car_deliv_pct),

            comp_justeat_bool = mean(comp_justeat_bool),
            comp_uber_bool = mean(comp_uber_bool),
            comp_justeat_rx_supply = round(mean(comp_justeat_rx_supply)),
            comp_uber_rx_supply = round(mean(comp_uber_rx_supply)),
            # num_rx_g6 = mean(num_rx_g6),
            # num_rx_enterprise = mean(num_rx_enterprise),
            # num_rx_local = mean(num_rx_local),
            num_rx_g25_pct = mean(num_rx_g25_pct),
            num_rx_enterprise_pct = mean(num_rx_enterprise_pct),
            #num_rx_local_pct = mean(num_rx_local_pct),
            #num_rx_nosegment_pct = mean(num_rx_nosegment_pct),
            # aof = mean(aof),
            # aov = mean(aov),
            #plus_orders_pct = mean(plus_orders_pct),
            # convenience_orders_pct = mean(convenience_orders_pct),
    
            # bicycle_deliv_pct = mean (bicycle_deliv_pct),
            # scooter_deliv_pct = mean(scooter_deliv_pct),
            # null_deliv_pct = mean(null_deliv_pct),
            # num_alcohol_rx_pct   = mean(num_alcohol_rx_pct),
            # num_american_rx_pct = mean(num_american_rx_pct),
            # num_british_rx_pct = mean(num_british_rx_pct ),
            # num_chinese_rx_pct = mean(num_chinese_rx_pct), 
            # num_dessert_rx_pct = mean(num_dessert_rx_pct),
            # num_french_rx_pct = mean(num_french_rx_pct),
            # num_greek_rx_pct = mean(num_greek_rx_pct),
            # num_indian_rx_pct  = mean(num_indian_rx_pct ),
            # num_italian_rx_pct = mean(num_italian_rx_pct),
            # num_japanese_rx_pct = mean(num_japanese_rx_pct), 
            # num_korean_rx_pct  = mean(num_korean_rx_pct ),
            # num_lebanese_rx_pct  = mean(num_lebanese_rx_pct), 
            # num_mexican_rx_pct = mean(num_mexican_rx_pct ),
            # num_spanish_rx_pct = mean(num_spanish_rx_pct),
            # num_thai_rx_pct = mean(num_thai_rx_pct),
            # num_turkish_rx_pct = mean(num_turkish_rx_pct),
            # num_vietnamese_rx_pct = mean(num_vietnamese_rx_pct),
            # num_missing_rx_pct = mean(num_missing_rx_pct)
            )
  
nums <- vapply(df_summ, is.numeric, FUN.VALUE = logical(1))
df_summ[,nums] <- round(df_summ[,nums], digits = 2)

df_compare <- as.data.frame(t(df_summ))
names(df_compare) <- c('Cluster 1', 'Cluster 2', 'Cluster 3', 'Cluster 4', 'Cluster 5')
df_compare <- df_compare[-1,]

df_compare

NA

# 
# cluster_df <- do.call(cbind, cluster_list)
# 
# cluster_df <- cluster_df[,c(1,2,6)]
# names(cluster_df)[1] <- "zone_code"
# cluster_df$zone_change <- with(cluster_df, ifelse(cust_5.cluster == nocust_5.cluster, "N", "Y"))
# 
# cluster_df_Y <- cluster_df[cluster_df$zone_change == 'Y',]
# 
# 
# change1 <- cluster_df_Y$zone_code[cluster_df_Y$cust_5.cluster == 3 & cluster_df_Y$nocust_5.cluster == 4]
# change2 <- cluster_df_Y$zone_code[cluster_df_Y$cust_5.cluster == 4 & cluster_df_Y$nocust_5.cluster == 3]
# change3 <- cluster_df_Y$zone_code[cluster_df_Y$cust_5.cluster == 4 & cluster_df_Y$nocust_5.cluster == 1]
# change4 <- cluster_df_Y$zone_code[cluster_df_Y$cust_5.cluster == 3 & cluster_df_Y$nocust_5.cluster == 5]

There was some debate around whether customer-related factors like AOF and number of new customers should be included in the clustering algorithm. We removed customer-related factors to see how the clusters changed. Only 7.69% of clusters changed.

Bootstrapping Kmeans

We need to understand how well our clustering algorithm has performed - do our 5 clusters represent actual structure in the data, or are they products of our algorithm? Kmeans is an unsupervised machine learning algorithm, so there is no “right” answer to compare our outcomes to. Furthermore, the clusters specified can vary depending on where the initial centroid (centre of the cluster) is set.

We check whether our clusters represent a true structure by seeing if they remain stable under plausible variations in the dataset. We check the stability of our clusters the following way:

1. Bootstrap a sample of our feature dataset
   • Sample with replacement
     
     

2. Implement the k-means algorithm on the bootstrapped data

3. Compare the similarity of the clusters in (2) to our initial clustering output
   • Look at pairwise combinations of zones, and compute how many remain in the same cluster
     
     

4. Iterate over steps 1 - 3, 500 times to see the cumulative proportion of zone pairs that dissolve over 500 iterations

For example, imagine that our original clustering results in a cluster:
\(A = \{EDN, CAM, CHE, TON\}\)

We obtain 10 possible pairwise combinations: \(EDN-EDN, EDN-CAM, EDN-CHE, EDN-TON, CAM-CAM\) \(CAM-CHE, CAM-TON, CHE-CHE, CHE-TON, TON-TON\)
Note that we account for repeats of the same zone, and these combinations might occure in the bootstrapped data.

Our first bootstrapped implementation throws up a cluster:
\(B =\{CAM, CAM, CHE, WCV\}\)

Note that some repetition occurs due to us sampling with replacement (CAM), which also causes EDN and TON to not be picked in the random sample. This cluster also has a new zone, WCV, which was not assigned to that cluster originally.
We obtain 4 possible pairwise combinations: \(CAM-CAM, CAM-CHE, CAM-WCV, CHE-WCV\). 2 of these 4 pairs occur in our original cluster, so we calculate our similarity index as \(Similarity = 2/4 = 0.5\)

We initially used the Jaccard index to evaluate the stability of our clusters, as that is the most popular method. However, we ultimately realised that this was not applicable to our business problem.

In the example above, we would calculate the Jaccard index as \(J(A,B) = \{A ∩ B\}/\{A U B\} = \{CAM,CHE\} / \{EDN, CAM, CHE, TON, WCV\} = 0.4\)

This would result in a low Jaccard Index, simply because zones EDN and TON were not sampled in the bootstrapped data. Sampling with replacement meant that our bootstrapped clusters would almost always have less zones than our original clusters. So, the Jaccard index would unfairly penalize us for that.

start <- Sys.time()

cluster_list <- list()
cluster_df_list <- list()
pairwise_list_boot <- list()


sample_size <- nrow(df_vars2)
n_iter <- 1:200
n_clust <- 1:k
  
for (i in n_iter) {

  rows_to_sample <- sample(c(1:sample_size), replace = T)
  df_boot <- df_vars2[rows_to_sample,]
  
  #df_boot <- df_vars2
  
      if(i == 1) {
        op.kmeans.boot  <- op.kmeans #original clusters
      } else {
      set.seed(sample(c(1:1000),1))
      op.kmeans.boot <- kmeans(df_boot[,2:ncol(df_boot)], k, iter.max = 10)
      }
  
  
  zone_code_names <- rownames(as.data.frame(op.kmeans.boot$cluster))
  zone_code_names <- sub(".[1-9]+", "", zone_code_names) #to prevent re-naming of repeated zones
  
  clusters <- data.frame(zone_code = zone_code_names,
                         cluster = op.kmeans.boot$cluster,
                         iteration = i)
  clusters$zone_code <- as.character(clusters$zone_code)
  
  pairwise_list_boot.sub <- list()
  
  for (m in n_clust) {
    sub.clust <- clusters$zone_code[clusters$cluster == m]
    for (n in 1:length(sub.clust)) {
      for (p in 1:length(sub.clust)) {
      pair <- paste(sort(c(sub.clust[n], sub.clust[p])), collapse = "-")
      pairwise_list_boot.sub[[pair]] <- pair
    }
    }

  }

      # for (m in 1:nrow(clusters)) {
      #   for( n in 1:nrow(clusters)) {
      #     if (clusters[m,2] == clusters[n,2]) {
      #       pair <- paste(sort(c(clusters[m,1], clusters[n,1])), collapse = "-")
      #       pairwise_list_boot.sub[[pair]] <- pair
      #     }
      #   }
      # }
  
  pairwise_list_boot[[length(pairwise_list_boot)+ 1]] <- pairwise_list_boot.sub
  
  for (n in n_clust) {
    cluster_list[[paste0(i, "-", n)]] <- as.vector(clusters$zone_code[clusters$cluster == n])
  }
  cluster_df_list[[i]] <- clusters

  paste0("Completed: Iteration ",i)
}

end <- Sys.time()
  
time.taken <- end - start  
  

df_final_boot <- do.call(rbind, cluster_df_list)

numer_list <- list()
denom_list <- list()

for (i in 2:length(pairwise_list_boot)) {
  orig <- pairwise_list_boot[[1]]
  test <- pairwise_list_boot[[i]]
  numer <- sum(test %in% orig)
  denom <- length(test)
  numer_list [[length(numer_list )+1]] <- numer
  denom_list [[length(denom_list )+1]] <- denom
  }

numer_df <- t(data.frame(numer_list))
denom_df <- t(data.frame(denom_list))
iter_df <- data.frame(t = c(1:length(numer_df)),
                      orig_pairs = numer_df,
                      total_pairs = denom_df)

iter_df$cumu_orig_pairs <- cumsum(iter_df$orig_pairs)
iter_df$cumu_total_pairs <- cumsum(iter_df$total_pairs)
iter_df$mismatch <- with(iter_df, 1- (cumu_orig_pairs/cumu_total_pairs))

iter_df_kmeans <- iter_df

We see that the stability of our clusters is quite high. Over time, around 13.08% of our zone-pairs end up segmented into different clusters. This means that 86.92% of our zone-pairs remain in the same cluster. This is great.

plot_ly( x = iter_df_kmeans$t, y = iter_df_kmeans$mismatch, mode = 'lines', type = 'scatter') %>%
  layout(title = "% of incorrect pairwise matches over time",
         xaxis = list(title = "Number of Iterations"), 
         yaxis = list(title = "% Incorrect Matches", range = c(0,0.4)))

# plot_ly( y = iter_comp$pct_dissolved, mode = 'lines', type = 'scatter') %>%
#   layout(title = "% of dissolved clusters over time",
#          xaxis = list(title = "Number of Iterations"),
#          yaxis = list(title = "% Dissolved Clusters"))

Final Outputs

Appendix

Clustering using Hierarchical Clustering

We ultimately decided against hierarchical clustering for 2 reasons:

• 1. Hierarchical clustering does not work well with large datasets
• 2. The final cluster outputs were not as intuitive as the kmeans clusters

We also use hierarchical clustering to segment our zones, to see how the results compare to to Kmeans. We use the same set of features used for Kmeans clustering, and start clustering zones based on how similar they are (calculated using Euclidean distance). We use Maximum/ Complete linkage clustering as this tends to produce more compact clusters.

set.seed(926)

df_vars3 <- df_vars2
df.dist <- dist(df_vars3[,2:ncol(df_vars3)], method = 'euclidean')
hc <- hclust(df.dist, method = 'ward.D2')

hcd <- as.dendrogram(hc)

plot(hc, labels = df_vars3$zone_code,hang = -1, cex = 0.4, xlab = "Zones")

op = par(bg = "#d7e0df")

plot(cut(hcd, h = 1.75)$upper, main = "upper tree",
     col.main = "#45ADA8", col.lab = "#7C8071", 
     col.axis = "#F38630", lwd = 3, lty = 3, sub = "")
# add axis
axis(side = 2, at = seq(0, 400, 100), col = "#F38630", labels = FALSE,
    lwd = 2)

# add text in margin
mtext(seq(0, 400, 100), side = 2, at = seq(0, 400, 100), line = 1,
    col = "#A38630", las = 2)

NA
NA
NA
NA

The full dendrogram is too large for discerning anything meaningful, so we display the upper part of the tree, where h > 2. The optimal cutoff height seems to be around 4 - this gets us 5 clusters. Cutting the dendrogram any higher would result in 3 clusters (too little).

Note: We did attempt kmeans with 6 clusters, but this ended up splitting Cluster 1 (Momentum - New Launch and Relaunched) into 2 segments. We settled on 5 clusters, as it did not make business sense to have 2 clusters catering to new zones.

sub_grp <- cutree(hc, k = 5)

df_final2 <-  cbind(df, cluster = sub_grp)
fviz_cluster(list(data = df_vars3[,2:ncol(df_vars3)], cluster = sub_grp), labelsize = 10)

df_summ2 <- df_final2 %>%
  group_by(cluster) %>%
  summarise(number_of_zones = n(),
            zone_tenure_days = round(mean(zone_tenure_days)),
            zone_population = round(mean(zone_population)),
            #area = mean(area),
            popn_density = round(mean(popn_density)),
            zones_in_city = round(mean(zones_in_city)),
            live_rx_total = round(mean(live_rx_total)),
            number_cuisines = mean(number_cuisines),
            #rx_per_thousand_popn = mean(rx_per_thousand_popn),
            has_editions_site = mean(has_editions_site),
            pct_mplus = mean(pct_mplus),
            avg_price_category = mean(avg_price_category),
            car_deliv_pct = mean(car_deliv_pct),

            comp_justeat_bool = mean(comp_justeat_bool),
            comp_uber_bool = mean(comp_uber_bool),
            comp_justeat_rx_supply = round(mean(comp_justeat_rx_supply)),
            comp_uber_rx_supply = round(mean(comp_uber_rx_supply)),
            # num_rx_g6 = mean(num_rx_g6),
            # num_rx_enterprise = mean(num_rx_enterprise),
            # num_rx_local = mean(num_rx_local),
            num_rx_g25_pct = mean(num_rx_g25_pct),
            num_rx_enterprise_pct = mean(num_rx_enterprise_pct),
            #num_rx_local_pct = mean(num_rx_local_pct),
            #num_rx_nosegment_pct = mean(num_rx_nosegment_pct),
            # aof = mean(aof),
            # aov = mean(aov),
            #plus_orders_pct = mean(plus_orders_pct),
            # convenience_orders_pct = mean(convenience_orders_pct),
    
            # bicycle_deliv_pct = mean (bicycle_deliv_pct),
            # scooter_deliv_pct = mean(scooter_deliv_pct),
            # null_deliv_pct = mean(null_deliv_pct),
            # num_alcohol_rx_pct   = mean(num_alcohol_rx_pct),
            # num_american_rx_pct = mean(num_american_rx_pct),
            # num_british_rx_pct = mean(num_british_rx_pct ),
            # num_chinese_rx_pct = mean(num_chinese_rx_pct), 
            # num_dessert_rx_pct = mean(num_dessert_rx_pct),
            # num_french_rx_pct = mean(num_french_rx_pct),
            # num_greek_rx_pct = mean(num_greek_rx_pct),
            # num_indian_rx_pct  = mean(num_indian_rx_pct ),
            # num_italian_rx_pct = mean(num_italian_rx_pct),
            # num_japanese_rx_pct = mean(num_japanese_rx_pct), 
            # num_korean_rx_pct  = mean(num_korean_rx_pct ),
            # num_lebanese_rx_pct  = mean(num_lebanese_rx_pct), 
            # num_mexican_rx_pct = mean(num_mexican_rx_pct ),
            # num_spanish_rx_pct = mean(num_spanish_rx_pct),
            # num_thai_rx_pct = mean(num_thai_rx_pct),
            # num_turkish_rx_pct = mean(num_turkish_rx_pct),
            # num_vietnamese_rx_pct = mean(num_vietnamese_rx_pct),
            # num_missing_rx_pct = mean(num_missing_rx_pct)
            )
  
nums <- vapply(df_summ2, is.numeric, FUN.VALUE = logical(1))
df_summ2[,nums] <- round(df_summ2[,nums], digits = 2)

test <- as.data.frame(t(df_summ2))
names(test) <- c('Cluster 1', 'Cluster 2', 'Cluster 3', 'Cluster 4', 'Cluster 5')
test <- test[-1,]

test

NA
NA

We now bootstrap our dataset to see how stable our clusters remain over time.

start <- Sys.time()

cluster_list <- list()
cluster_df_list <- list()
pairwise_list_boot <- list()

k_hc <- 5
sample_size <- nrow(df_vars3)
n_iter <- 1:250
n_clust <- 1:k_hc
  
for (i in n_iter) {

  rows_to_sample <- sample(c(1:sample_size), replace = T)
  df_boot <- df_vars3[rows_to_sample,]
  
  #df_boot <- df_vars2
  
    if(i == 1) {
      df.dist.boot <- df.dist
      hc.boot <- hc

    } else {
      set.seed(sample(c(1:1000),1))
      df.dist.boot <- dist(df_boot[,2:ncol(df_boot)], method = 'euclidean')
      hc.boot  <- hclust(df.dist.boot, method = 'ward.D2')

  }
  
  sub_grp.boot <- cutree(hc.boot, k = k_hc)
  
  zone_code_names <- hc.boot$labels
  zone_code_names <- sub(".[1-9]+", "", zone_code_names) #to prevent re-naming of repeated zones
  
  clusters <- data.frame(zone_code = zone_code_names,
                         cluster = sub_grp.boot,
                         iteration = i)
  clusters$zone_code <- as.character(clusters$zone_code)

  
  pairwise_list_boot.sub <- list()
  
  for (m in n_clust) {
    sub.clust <- clusters$zone_code[clusters$cluster == m]
    for (n in 1:length(sub.clust)) {
      for (p in 1:length(sub.clust)) {
      pair <- paste(sort(c(sub.clust[n], sub.clust[p])), collapse = "-")
      pairwise_list_boot.sub[[pair]] <- pair
    }
    }

  }

      # for (m in 1:nrow(clusters)) {
      #   for( n in 1:nrow(clusters)) {
      #     if (clusters[m,2] == clusters[n,2]) {
      #       pair <- paste(sort(c(clusters[m,1], clusters[n,1])), collapse = "-")
      #       pairwise_list_boot.sub[[pair]] <- pair
      #     }
      #   }
      # }
  
  pairwise_list_boot[[length(pairwise_list_boot)+ 1]] <- pairwise_list_boot.sub
  
  for (n in n_clust) {
    cluster_list[[paste0(i, "-", n)]] <- as.vector(clusters$zone_code[clusters$cluster == n])
  }
  cluster_df_list[[i]] <- clusters

  paste0("Completed: Iteration ",i)
}

end <- Sys.time()
  
time.taken <- end - start  
  

df_final_boot <- do.call(rbind, cluster_df_list)

numer_list <- list()
denom_list <- list()

for (i in 2:length(pairwise_list_boot)) {
  orig <- pairwise_list_boot[[1]]
  test <- pairwise_list_boot[[i]]
  numer <- sum(test %in% orig)
  denom <- length(test)
  numer_list [[length(numer_list )+1]] <- numer
  denom_list [[length(denom_list )+1]] <- denom
  }

numer_df <- t(data.frame(numer_list))
denom_df <- t(data.frame(denom_list))
iter_df <- data.frame(t = c(1:length(numer_df)),
                      orig_pairs = numer_df,
                      total_pairs = denom_df)

iter_df$cumu_orig_pairs <- cumsum(iter_df$orig_pairs)
iter_df$cumu_total_pairs <- cumsum(iter_df$total_pairs)
iter_df$mismatch <- with(iter_df, 1- (cumu_orig_pairs/cumu_total_pairs))

iter_df_hclust <- iter_df

We see that our hierarchical clusters are also highly stable, and actually marginally more stable than to our kmeans clusters. Over time, around 11.7% of our zone-pairs end up segmented into different clusters. This means that 88.3% of our zone-pairs remain in the same cluster. While hierarchical clustering gave us higher cluster stability, we decided to go with kmeans for the following reasons: * The 2 methodologies led to basically the same clusters! * The difference in performance was marginal (1.38 pp) in terms of cluster stability *Hierarchical clustering is more computationally expensive

plot_ly( x = iter_df_hclust$t, y = iter_df_hclust$mismatch, mode = 'lines', type = 'scatter') %>%
  layout(title = "% of incorrect pairwise matches over time",
         xaxis = list(title = "Number of Iterations"), 
         yaxis = list(title = "% Incorrect Matches", range = c(0,0.3)))

# cb1 <- clusterboot(df_vars2[2:ncol(df_vars2)], B = 100, bootmethod = 'boot', bscompare = T,
#                    multipleboot = T, clustermethod = kmeansCBI, dissolution = 0.5, krange = 5,
#                    seed = 18)
# 
# df_cb <- cb1$result$result
# groups<-cb1$result$partition
# cb1$bootmean
# cb1$bootbrd


# 
# iter_l <- data.frame(cluster_first = integer(),
#                    cluster_second = integer(),
#                    J.dist = numeric())
# 
# l <- 1:length(cluster_list)
# 
# for (i in n_clust) {
#   x <- cluster_list[[i]]
#   
#   for(c in (l[!(l %in% i)]) ) {
#     y <- cluster_list[[c]]
# 
#     t <- data.frame(cluster_first = names(cluster_list)[[i]], 
#                    cluster_second = names(cluster_list)[[c]],
#                    J.dist = length(unique(x[x %in% y]))/length(unique(c(x, y)))       )
#     
#      
#     iter_l <- rbind(iter_l,t)
#     
#   }
#   
# }
# 
# 
# 
# iter_l$iter_x <- as.numeric(sub("-.+", "", iter_l$cluster_first))
# iter_l$clus_x <- as.numeric(sub(".+-", "", iter_l$cluster_first))
# 
# iter_l$iter_y<- as.numeric(sub("-.+", "", iter_l$cluster_second))
# iter_l$clus_y <- as.numeric(sub(".+-", "", iter_l$cluster_second))
# 
# iter_summ <- iter_l %>%
#   group_by(iter_x, clus_x, iter_y) %>%
#   mutate(max_Jdist = max(J.dist)) 
# 
# iter_summ <- iter_summ %>%
#   filter(J.dist == max_Jdist) %>%
#   filter(iter_x != iter_y) %>%
#   mutate(cluster_dissolved = ifelse(J.dist < 0.5, 1, 0))
# 
# iter_summ <- iter_summ[order(iter_summ$iter_x, iter_summ$iter_y, iter_summ$clus_x),]
#          
# iter_summ$cum_sum_cluster_dissolved <- cumsum(iter_summ$cluster_dissolved)
# iter_summ$num <- 1
# iter_summ$total_clusters <- cumsum(iter_summ$num)
# 
# iter_comp <- iter_summ[iter_summ$clus_x == k,] #ie obtain the cumulative number of dissolved clusters, and total clusters, at the end of each iteration
# 
#   
# iter_comp$pct_dissolved <- iter_comp$cum_sum_cluster_dissolved/iter_comp$total_clusters